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Abstract. - The electronic structure of the monoclinic structure of Fe304 is studied using 
both the local density approximation (LDA) and the LDA-|-[/. The LDA gives only a small 
charge disproportionation, thus excluding that the structural distortion should be sufficient to 
give a charge order. The LDA-I-J7 results in a charge disproportion along the c-axis in good 
agreement with the experiment. We also show how the effective U can be calculated within 
the augmented plane wave methods. 



Introduction. - Magnetite has received a lot of attention both for fundamental and 
technological reasons. Above the Verwey transition it is a half-metal and has the highest 
known of 860 K. The crystal structure is the inverse spinel structure with the formal 
chemical composition Fe^'''[Fe^+Fe'^+]B04~. The two octahedrally coordinated B positions 
are symmetry equivalent and order antiferromagnetically with the tetrahedrally coordinated A 
site in the cubic Fd3m spacegroup. When cooled to the Verwey transition temperature, which 
lies around 122-125 K depending on sample, the conductivity of magnetite drops abruptly by 
two orders of magnitude. Originally the structure below the Verwey transition transition 
was thought to be have iron cations at the B sites order as Fe^"*" and Fe'^+ along the (001) 
planes. This turned out to be too simple a model and single crystal diffraction studies showed 
that the low temperature structure is monoclinic with a x V^a x 2a supercell and a 
Cc space group symmetry. A recent NMR study resolved 8 tetrahedral and 16 octahedral 
environments thereby confirming the Cc space group. [1] However, the diffracted supercell 
peaks are extremely weak and even a recent synchrotron diffraction study could only resolve 
three supercell peaks. [2] The structure has therefore only been refined in a a /%/2 X a/V2x 2a 
monoclinic subcell with an additional orthorhombic Pmca pseudo symmetry constraint (see 
refs. [2,3] for a detailed description of the structure). 

It is generally agreed that the charge order model with alternating Fe^~'"/Fe'^+ layers below 
the Verwey transition is too simple but the full charge order is still not fully understood. A 
Mossbauer study fitted the spectrum with five components. [4] One corresponding to Fe^"^ and 
four to Fe^"^ and Fe'^"^ on two non-equivalent octahedral sites. [4] A resonant X-ray diffraction 
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study reported that if a charge disproportionation takes place it should be below the sensitivity 
limit of 25% of the experiment. [5] Recently it was noticed that the structural distortions lead 
to significantly different mean Fe— O distances for the octahedrally coordinated iron sites. [2] 
When these distances are used to calculate bond valence sums a charge modulation along 
[001] is found and the octahedral iron sites split into two groups with a charge disproportion 
of about 0.2 electron. Furthermore a [OO5] second charge modulation is found, which leads to 
the doubling of the c-axis compared to the cubic structure. [2] 

Theoretically the problem of charge order in Fe304 has been difficult to attack. First 
of all the true structure is unknown and even the simplified a/\/2 x a/\/2 x 2a subcell is 
highly complex for a theoretical calculation. Furthermore, the local density approximation 
(LDA) is not generally applicable to highly correlated transition metal oxides, due to the 
spurious self interaction. F.inst. if the calculations are done in the originally proposed Verwey 
model without structural distortions, an LDA calculation converges to a metallic state with 
the octahedral iron sites in an Fe^'^"*" oxidation state. [6] In the LDA+[/ method an orbital 
dependent field is introduced which can be shown to give a correction for the self interaction. [7] 
Consequently two studies have shown that the LDA+J7 method leads to a charge ordering. [6,8] 
However, both these studies used the original Verwey model of the structure thereby omitting 
both the structural distortions and the additional [00^] charge modulation. [2] Till date only 
one calculation has been done on the a/\/2 x a/\/2 x 2a structure. [9] Both a pure LDA 
and three self interaction corrected (SIC)-LDA calculations were performed and the authors 
reached two surprising conclusions. [9] First of all it was found that the scenario where five 
d-electrons move in the SIC-LDA potential on both the A and all B sites (which can be 
interpreted as all the octahedral iron sites being in the Fe^"*" oxidation state) was the most 
stable. [9] Secondly the lowest energy state shows only a small charge disproportion of 0.1 
electron. This is of similar magnitude to the charge disproportion found using pure LDA, 
leading to the conclusion that the charge disproportion is structural of origin. [9] One possible 
objection to these results could be that the minimal basis set and atomic sphere approximation 
used for the potential [9] could make the total energies unreliable. A second objection could 
be that in the SIC-LDA method one must choose a set of localized states for which the SIC 
is applied. The actual ground state will thus only be found if it is among those tested. 

We have therefore performed LDA+J7 calculations using the linearized augmented plane 
wave(LAPW) method. Furthermore we have tested a method for deriving i7-values using the 
LAPW method. We do indeed find an electronic configuration that was not considered in 
the SIC-LDA study [9], but which is in good agreement with the charge disproportion that 
was proposed based on the experimental bond distances. [2] The paper is arranged into two 
sections. First we describe the LDA-[/ method and how we derive [/-values within the LAPW 
method. Secondly the calculated electronic structure is described. 

Fixing the U within the LAPW method. - The LDA+J/ method essentially consists of 
identifying a set of atomic like orbitals which are treated in a non-LDA manner. [10] Based on 
the lessons from Hubbard model studies these orbitals are treated with an orbital dependent 
potential with an associated on-site Coulomb and exchange interactions, U and J. To avoid 
double counting in the non-spherical part of potential we use C/e// = U — J instead of U [11] 
and omit the multipolar terms proportional to J in the added LDA-l-17 potential. 



where h is the orbital occupation matrix. This form is both rotationally invariant and is well 
suited for full potential calculations. 




(1) 
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In LDA the electron-electron interactions have already been taken into account in a mean 
field way. One must therefore identify the parts that occur twice and apply a double counting 
correction (DCC). Several different versions of the DCC exist. [7,10,14,15] In the present work 
we will use what has been called the fully localized limit (FLL). [7] The FLL-DCC has a clear 
interpretation and has also been shown to give results in better agreement with experiment 
than other DCC. [16]. The FLL-DCC is constructed so that the atomic like limit to the total 
energy is satisfied, which leads to the following expression [7] 

E^S^ = -i^nin - 1) _ ^ ^ - 1)) . -^^.(^1 + l)n^ (2) 

where ria = Tr{n'^)/ {21 -\- \). The orbital dependent potentials entering the Kohn-Sham 
equation that arise from the E""^^ — pDCC correction to the total energy, Eq. Q and Eq. (O 
is then 

^yu.^ = qJ^'^^^ - -{U - J){n^ \r) (3) 

which shows that the FLL version of the orbital potential will stabilize an orbital that is more 
than half occupied and destabilize an orbital that is less than half occupied. 

The meaning of the U parameter was discussed by Anisimov and Gunnarsson, [12] who 
defined it as the cost in Coulomb energy by placing two electrons on the same site. In an atom 
the U corresponds to F° of the unscreened Slater integrals. [12] F° should thus both increase 
with increased ionicity and as the d-wave function is contracted across the 3d transition series. 
This is illustrated in Fig. ^where we show the calculated atomic Slater integrals for chemically 
relevant 3d ions. 

Due to screening the effective U in solids is much smaller than F'^ for atoms. To calculate 
the effective U Anisimov and Gunnarsson, [12] constructed a supercell and set the hopping 
integrals connecting the 3d orbital of one atom with all other orbitals to zero. The number of 
electrons in this non-hybridizing d-shell was varied and -Fg// then calculated from 

F°ff ^e^d^{{n + l)/2,n/2) - e^d^({n + l)/2, n/2 - 1) 
- SFiin + l)/2, n/2) + £j.((n + l)/2, n/2 - 1) 

where £3^1 is the spin-up 3d eigenvalue. Using the method of Anisimov and Gunnarsson, 
/7-values have been calculated for the di- and trivalent configurations of the 3d elements in 
La perovskites. [13] As expected the trends were the same as observed for atoms. Fig. but 
in smaller magnitude. The trends being, i), an almost linear relation between atomic number 
and calculated U and ,ii), a constant shift between the M^"*" ions, ranging from approximately 
6.5 eV (Titanium) to 8.5 eV (Copper), and the ions, ranging from 8 eV to 10 eV. [13] 

The original LDA-I-J7 method [10] was based on the linearized muffin tin orbitals basis set, 
where the individual orbital and hopping terms can be identified. This is not possible within 
the LAPW method, so the method of Anisimov and Gunnarsson [12] can not be directly 
applied. Instead the hybridization can be removed putting the d-states into the core or by 
performing a two- window calculation. To check this procedure we have performed calculations 
on the well characterized NiO. [11] Two calculations were performed on 2 x 2 x 2 supercells each 
with one impurity site with the d-configuration forced to be as in Eq. Q . The d-character of 
the APWs at the impurity sites was eliminated by placing the d-linearization energy at a very 
high value. Using Eq. Q we hereby got a value of Fgj^=5.96 eV. The question is how this value 
should be used in an LDA-|-t/ calculation? When using the spherically averaged form of E""^^ , 
Eq. J simply rescales the orbital term, Eq. (O and Eq. 10). Furthermore Eq. ^ shows 
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Ti V Cr Mn Fe Co Ni Cu Zn Ti V Cr Mn Fe Co Ni Cu Zn 



Fig. 1 - Onsite parameters calculated from the atomic Slater integrals. □ atomic values, A M^"*" 
ions and o M''"'" ions. The atomic values for Cr and Cu were calculated with 5 and 10 d-electrons 
respectively. 



that an occupied and an unoccupied orbital will be split hy U — J. Following the interpretation 
of Anisimov and Gunnarsson [12], this suggest that U — J — F^ff- As the screening of 
and in solids appears to be small, J can to a good approximation be calculated from the 
atomic values. [12] Using a J = 1.36 eV, Fig.QJ we arrive at U — 7.32 eV in good agreement 
with earlier values. [11] Furthermore an effective onsite term of F^^j—b.^Q eV, leading to 
a corresponding splitting of occupied and unoccupied bands, is in excellent agreement with 
experiment [11] 

We then applied the method to FeO and Fe203 and obtained — 5-73 eV and 7.33 
eV for Fe^+ and Fe^+ respectively. Finally calculations were performed on magnetite using 
the 2x2x2 supercell of the cubic high temperature structure. We tried placing both Fe^'*' 
(corresponding to one calculation with 3.5 spin up and 3 spin down electrons at the impurity 
site and one calculation with 3.5 spin up and 2 spin down) and Fe^"*" at the tetrahedrally 
coordinated A site. Thereby values of =4.79 eV and 6.33 eV were obtained. For the 
octahedrally coordinated B site values of F0^^(Fe|+)=5.03 eV, F0^^(Fe^^+)=6.21 eV and 
Fgjj(Fe^"^)=7.38 eV were calculated, showing good internal consistency with the FeO and 
Fe2 03 values. Using the J parameters derived from the atomic values, Fig.^ we then obtain 
[/(Fe^+)=7.69 eV, t7(Fe|+)=6.2 eV and f/(Fe|j+)=8.73 eV for magnetite. Our calculated F^^^ 
seem somewhat larger than the U and J values that were earlier calculated for magnetite. 
[6] However it is not clear what oxidation state was used to calculate these values [6] and 
furthermore it should be noted that the values for the octahedrally coordinated B site are in 
very good agreement with what was obtained for the octahedrally coordinated Fe ions in the 
perovskites. [13] 

Computational details. - With the U values thus fixed we performed calculations on the 
low-symmetry structure [2] using the L/APW+lo method [17], as implemented in the WIEN2k 
code. [18] Sphere sizes of 2.0 a.u. and 1.5 a.u. were used for the Fe and O sites respectively. 
The plane-wave (PW) cut-off was defined by min(i?MT)max(fc„)=5.7 corresponding to ap- 
proximately 3700 PW. The Brillouin-zone (BZ) was sampled on a tetrahedral mesh with 100 
k points (18 in the irreducible BZ). On the tetrahedrally coordinated A iron sites the calcu- 
lated F°ff ~ J ^ 6.33 eV was used. The calculated F^^^ = 6.21 eV for Fe^^+ was used 
on all octahedrally coordinated B iron sites so that no charge order was forced. 

When discussing charge order it is necessary to define some kind of "atomic" charge. In the 
APW method these partial charges are most commonly the integrated charge within a given 
sphere. They can be named q*, where t refers to a given atomic sphere and I to the angular 
quantum number. The values have the problem that they depend on the size of the atomic 
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2.30 
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Table I - Fractional charges, valencies and Spin (nig) magnetic moments (in fis) calculated as 
explained in the computational section. The primed values have been scaled by a factor 1.06 (see the 
computational section for explanation). The sites are labeled as in Ref. [2]. The Bl and B2 sites are 
both an average over two symmetry-nonequivalent sites. These sites would be symmetry equivalent in 
an orthorhombic unit cell and as expected (the mono-clinic distortion is very small, (3 = 90.2363 [2]^ 
the values calculated for these sites are very similar and therefore only one averaged value is reported. 



sphere, thus depend on a computational parameter. They are still useful for comparing atomic 
charges within the same calculation, if the sphere sizes are kept equal. They can however not 
directly be compared with the results of other calculations/experiments. One way to make 
the charges more transferable is to estimate the amount of charge leaking out of the sphere 
and use this for rescaling. To do the rescaling we have calculated the radial function of free 
Fe^+ and Fe"^+ ions. By integrating the d-functions up to 2 a.u. (the sphere size used for iron 
in the present study) we get an average scaling of 1.06. In the following two valencies will be 
defined for the iron sites = 8 — (the Qq and g* values are less than 0.05 electrons) and 
V*' = 8 — 1.06 X As will be shown in the next section the scaling seems to be reasonable 
as it brings fully localized half shells up to 5 electrons. 

Results and discussion. - Table Ogives the fractional charges for Fe304. It is seen that 
LDA gives only a very small charge disproportionation, thus excluding that the structural 
distortion should be sufficient to give a charge order. In contrast the LDA-I-J7 calculation 
does give a charge disproportionation, TableQlof approximately 0.2 electrons between the B2 
and B3 (high oxidation state) and the Bl and BA (low oxidation state) sites, thus confirming 
the additional [00^] charge modulation. [2] Tabled also reports the rescaled partial charges. 
It is seen that the scaling brings the charge in the majority spin d-states of the high oxidation 
state iron sites close to five. This is expected and lends credibility to our scaling procedure. 
The scaled LDA+U results in a charge disproportion of 0.2—0.3 electrons and the valencies 
are found to be smaller than what could be expected from a naive ionic picture. 

Figure El shows the calculated density of states (DOS) for Fe304 in the charge ordered 
phase. For the LDA calculation the oxygen band is found between approximately -4 eV and 
-8 eV. Above these bands one finds the bands of Fe-Sd origin. Furthermore a substantial 
hybridization between the iron sites and the oxygen p-states is seen. When the orbital field is 
introduced the fully occupied minority spin bands on the A-iron sites are shifted downwards 
below the oxygen p-shell. The same thing happens for the B2 and B3 iron sites. This was 
expected from Table 13 where it is seen that the B2 and B3 have a valency similar to the A 
sites and that the majority spin bands seem to be fully occupied. FigureElalso shows that the 
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n(E] [#states/(u.c. spin eV)] n(e) [#states/(u.c. spin eV]] n(K) [#states/(u.c. spin eV)] 

Fig. 2 - Spin polarized density of states calculated using LDA (top left) and FLL (top right). The 
lower plots illustrate the partial DOS. Lower left: the A (tetrahedral) sites. Lower middle: 51+54 
(octahedral with long B-O bonds) sites. Lower right: B2+B3 (octahedral with short B-O bonds) 
sites. 

picture is very different for the Bl and sites. These sites are the ones with the long Fe— O 
distances [2] and low oxidation states, Table |l| It is seen that there is a strong hybridization 
between the majority spin Fe-3d and the 0-2p states. Furthermore there is a non-hybridizing 
minority spin states just below the Fermi level which is not found for the Bl and BA states 
and is the cause of the charge order. The electronic structure of the low oxidation state irons 
sites can thus be compared to the electronic structure of FeO, [19] where the Fe-Sd states also 
show a strong hybridization with the 0-2p states and one sub-band in the minority spin is 
occupied to accommodate the additional electron. 

The experimental work [2] chose to renormalize their bond valence sums so that the average 
valency of the B sites was 2.5, which is close to the average unsealed valency. Tabled in the 
present study. So despite the derived charges being method dependent (see computational 
section) the charges should be comparable. Consequently the valencies quoted 2.39, 2.61, 2.59 
and 2.41 [2] (Bl, B2, B3 and BA site respectively) which are in extremely good agreement 
with the present. Table ^ This we believe strongly supports our electronic structure, but 
then the question arrises why the earlier computational study of the charge ordered magnetite 
structure [9] did not reach the same conclusion. 

The earlier computational work [9] reported atomic charges and spin magnetic moments 
for six different B sites in the monoclinic phase for four different types of calculations. One cal- 
culation (labeled SIC(l)) is performed where the simple Verwey charge order is implemented, 
so that the B sites alternate between five d-electrons moving in the SIC potential and six 
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d-electrons moving in a SIC potential. [9] In our labeling a simple Verwey charge order corre- 
sponds to the atoms Bl and B2 being in a low oxidation state. These sites are actually split 
into two sites each in the monoclinic symmetry (see comment to Table ^) and consequently 
the SIC study finds four sites (in this paper labeled Bl, B2, B3 and _B4 [9]) that have a low 
oxidation state, while two sites ( BS and BA in our notation, i?5 and B6 in their) have a 
high oxidation state. This immediately explains why these authors found the charge ordered 
state to be higher in energy than the state where all the B sites had five d-electrons moving 
in a SIC potential. The simple Verwey state considered earlier [9] does not correspond to the 
charge ordering that could be expected from the structural distortion, [2] which means that 
the B3 site is calculated to be in the high oxidation state despite it having longer distances. 
In other words the SIC calculation have omitted the second [00^] charge modulation found 
experimentally. [2] 

Conclusion. - By electronic structure calculations using the LDA+[/ method we have 
obtained a periodic charge disproportion along the c-axis in the monoclinic structure of Fe304. 
The proposed charge order is in good agreement with the one derived from experiment and 
we have explained the disagreement with an earlier SIC-LDA calculation. 

This work was supported by the project A1010214 of the Grant Agency of the AS CR and 
by the Carlsberg Foundation. 
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